program drop _all
program define DPDprog, rclass
version 14
syntax, beta0(real) betax(real) betay(real) numN(integer) numT(integer) numNT(integer) ///
        beffect(real)
drop _all
set obs `numN'
gen id = _n
// "ai" is the part of the error term that doesn't change over time for a given
// unit.
gen ai = rnormal()
expand `numT'
bysort id: gen t=_n
xtset id t
gen x = (`beffect'*ai + rnormal())/sqrt(1+`beffect'^2)
// "uit" is the part of the error term that varies across time
gen uit = rnormal()
// The total error term is the sum of "a" and "e"
gen error = ai + uit
gen y = `beta0'/(1-`betay')
replace y = `beta0' + `betax'*x + `betay'*L.y + error if t > 1


// This section calculates DFE estimates
regress y x L.y i.id if t > 20
return scalar DFELRP = _b[x]/(1-_b[L.y])
testnl _b[x]/(1-_b[L.y])=`betax'/(1-`betay')
return scalar pDFELRP = r(p)

//This section calculates AH estimates
xtivreg y x (L.y = L(2/3).y) if t > 20, fd
return scalar AHLRP = _b[D.x]/(1-_b[LD.y])
testnl _b[D.x]/(1-_b[LD.y])=`betax'/(1-`betay')
return scalar pAHLRP = r(p)

// This section calculates Difference GMM
xtabond y x if t > 20
return scalar DGMMLRP = _b[x]/(1-_b[L.y])
testnl _b[x]/(1-_b[L.y])=`betax'/(1-`betay')
return scalar pDGMMLRP = r(p)

// This section calculates System GMM
xtdpdsys  y x if t > 20
return scalar SGMMLRP = _b[x]/(1-_b[L.y])
testnl _b[x]/(1-_b[L.y])=`betax'/(1-`betay')
return scalar pSGMMLRP = r(p)

end
